
## PEASANT RESISTANCE AND ECONOMIC AFFLUENCE: LESSONS FROM PARAGUAY
## Appendix D: Extensions
## January 2024

## Code Authors:
## Germán Feierherd (gfeierherd@udesa.edu.ar) 
## Jorge Mangonnet (jorge.mangonnet@vanderbilt.edu)


#### Preamble ------------------------------------------------------------------

## Workspace setup
rm(list = ls())		         # clear all objects in memory
options(max.print = 5000)  # expand display
options(scipen = 999)      # disable sci notation
dev.off()                  # reload graphic device
cat("\014")                # clear console

## Load/install packages
if (!require("pacman")) install.packages("pacman")
pacman::p_load(
  readxl,
  interflex,
  marginaleffects,
  haven,
  readr,
  reshape2,
  magrittr,
  dplyr,
  plm,
  lfe,
  pcse,
  lmtest,
  stargazer,
  xtable,
  texreg,
  estimatr, 
  fixest,
  margins, 
  sjPlot, 
  multiwayvcov,
  #  rgdal,
  ggplot2,
  wesanderson, 
  fixest,
  did,
  sf)

## Set working directory
gf <- "~/Dropbox/Working_Papers/investigacion protesta Paraguay/"; main <- gf  # german's path
jm <- "~/Dropbox/Research/investigacion protesta Paraguay/"; main <- jm        # jorge's path
setwd(main)

## Functions
source(paste0(main, "4_code/functions_paraguay.R"))

## Load full dataset
load(paste0(main, "4_code/DataSetPY.Rda"))


#### Appendix D: Extensions ----------------------------------------------------

#### Table D.1: Base Results by Capital-Intensive Crop ----

# Model 1: Soybeans, w/o controls
m1 <- plm(ln_conflict ~ I(ln_soy_suit*ln_soy_price) + factor(year), index = c("codigo"), model = "within", effect = "individual", data = dta); summary(m1) 
se1 <- coeftest(m1, vcov = vcovHC(m1, type = "HC2", cluster = "group")); se1

# Model 2: Soybeans, w/ controls
m2 <- plm(ln_conflict ~ I(ln_soy_suit*ln_soy_price) + ln_pop + factor(year)*distance_hubs + factor(year)*untitled_share, index = c("codigo"), model = "within", effect = "individual", data = dta); summary(m2) 
se2 <- coeftest(m2, vcov = vcovHC(m2, type = "HC2", cluster = "group")); se2

# Model 3: Maize, w/o controls
m3 <- plm(ln_conflict ~ I(ln_mze_suit*ln_mze_price) + factor(year), index = c("codigo"), model = "within", effect = "individual", data = dta); summary(m3) 
se3 <- coeftest(m3, vcov = vcovHC(m3, type = "HC2", cluster = "group")); se3

# Model 4: Maize, w/ controls
m4 <- plm(ln_conflict ~ I(ln_mze_suit*ln_mze_price) + ln_pop + factor(year)*distance_hubs + factor(year)*untitled_share, index = c("codigo"), model = "within", effect = "individual", data = dta); summary(m4) 
se4 <- coeftest(m4, vcov = vcovHC(m4, type = "HC2", cluster = "group")); se4

# Model 5: Sugar, w/o controls
m5 <- plm(ln_conflict ~ I(ln_sug_suit*ln_sug_price) + factor(year), index = c("codigo"), model = "within", effect = "individual", data = dta); summary(m5) 
se5 <- coeftest(m5, vcov = vcovHC(m5, type = "HC2", cluster = "group")); se5

# Model 6: Sugar, w/ controls
m6 <- plm(ln_conflict ~ I(ln_sug_suit*ln_sug_price) + ln_pop + factor(year)*distance_hubs + factor(year)*untitled_share, index = c("codigo"), model = "within", effect = "individual", data = dta); summary(m6) 
se6 <- coeftest(m6, vcov = vcovHC(m6, type = "HC2", cluster = "group")); se6

# Create table 
stargazer(
  m1, m2, m3, m4, m5, m6,  
  se = list(se1[,2], se2[,2], se3[,2], se4[,2], se5[,2], se6[,2]), 
  p = list(se1[,4], se2[,4], se3[,4], se4[,4], se5[,4], se6[,4]),  
  covariate.labels = c("Soybean Prices$_{t-1}$ $\\times$ Soybean Suitability",
                       "Maize Prices$_{t-1}$ $\\times$ Maize Suitability",
                       "Sugar Prices$_{t-1}$ $\\times$ Sugar Suitability"),
  add.lines = list(c("","","","","","",""), 
                   c("Controls", rep(c("No", "Yes"), 3))),
  #c("Municipality FE", rep("Yes", 8)),
  #c("Year FE", rep("Yes", 8))),
  omit = c("factor", "Constant", "distance_hubs", "untitled_share", "ln_pop"),
  omit.stat = c("adj.rsq", "ser", "f"),
  dep.var.labels.include = FALSE,
  model.names = FALSE,
  df = FALSE, 
  multicolumn = TRUE, 
  star.cutoffs = c(0.1, 0.05, 0.01),
  star.char = c("*", "**", "***"),
  #float = TRUE, 
  #float.env = "sidewaystable",
  font.size = "footnotesize",
  column.sep.width = "-1.5pt",
  style = "ajps", 
  title = "Base Results by Capital-Intensive Crop", 
  label = "commodity",
  type = "latex",
  notes.append = FALSE,
  notes = c("\\parbox[t]{15.25cm}{\\textit{Note}: All models include municipality and year fixed effects. Standard errors clustered by municipality in parentheses.}",
            "$^{*}p<.1$, $^{**}p<.05$, $^{***}p<.01$"),
  out = paste0(main,"1_writing/2023/tables/table_commodity.tex")
)


#### Figure D.1: Marginal Effect of Agricultural Prices by Land Suitability for each Capital-Intensive Crop ----

# Based on model 1, 3 & 5 
m1 <- plm(ln_conflict ~ ln_soy_suit*ln_soy_price + factor(year), index = c("codigo"), model = "within", effect = "individual", data = dta); summary(m1) 
m3 <- plm(ln_conflict ~ ln_mze_suit*ln_mze_price + factor(year), index = c("codigo"), model = "within", effect = "individual", data = dta); summary(m3) 
m5 <- plm(ln_conflict ~ ln_sug_suit*ln_sug_price + factor(year), index = c("codigo"), model = "within", effect = "individual", data = dta); summary(m5) 

# Covariance matrix for cluster standard errors
covMat1 <- vcovHC(m1, type = "HC2", cluster = "group")
covMat3 <- vcovHC(m3, type = "HC2", cluster = "group")
covMat5 <- vcovHC(m5, type = "HC2", cluster = "group")

# Plots
pdf(file = paste0(main, "1_writing/2023/figures/inter_soybean.pdf"), width = 5, height = 5)
interaction_plot_continuous(m1, effect="ln_soy_price", moderator="ln_soy_suit", interaction="ln_soy_suit:ln_soy_price", varcov=covMat1, minimum=0.66, maximum=1.12, incr="default", num_points = 100, conf=.95, mean=FALSE, median=FALSE, alph=50, rugplot=F, histogram=T, title=NULL, xlabel="Soybean Suitability", ylabel="Marginal Effect of Soybean Prices")
dev.off()

pdf(file = paste0(main, "1_writing/2023/figures/inter_maize.pdf"), width = 5, height = 5)
interaction_plot_continuous(m3, effect="ln_mze_price", moderator="ln_mze_suit", interaction="ln_mze_suit:ln_mze_price", varcov=covMat3, minimum=1.40, maximum=2.4, incr="default", num_points = 100, conf=.95, mean=FALSE, median=FALSE, alph=50, rugplot=F, histogram=T, title=NULL, xlabel="Maize Suitability", ylabel="Marginal Effect of Maize Prices")
dev.off()

pdf(file = paste0(main, "1_writing/2023/figures/inter_sugar.pdf"), width = 5, height = 5)
interaction_plot_continuous(m5, effect="ln_sug_price", moderator="ln_sug_suit", interaction="ln_sug_suit:ln_sug_price", varcov=covMat5, minimum=-4.5, maximum=-4, incr="default", num_points = 100, conf=.95, mean=FALSE, median=FALSE, alph=50, rugplot=F, histogram=T, title=NULL, xlabel="Sugar Suitability", ylabel="Marginal Effect of Sugar Prices")
dev.off()

remove(m1, m2, m3, m4, m5, m6, se1, se2, se3, se4, se5, se6, covMat1, covMat3, covMat5)


#### Table D.2: Excluding Municipalities from the Chaco Region ----

# Model 1: Prices x Suitability, w/o controls
m1 <- plm(ln_conflict ~ I(ln_suitability*ln_price) + factor(year), index = c("codigo"), model = "within", effect = "individual", subset = region != "Chaco", data = dta); summary(m1) 
se1 <- coeftest(m1, vcov = vcovHC(m1, type = "HC2", cluster = "group")); se1

# Model 2: Prices x Suitability, w/ controls
m2 <- plm(ln_conflict ~ I(ln_suitability*ln_price) + ln_pop + factor(year)*distance_hubs + factor(year)*untitled_share, index = c("codigo"), model = "within", effect = "individual", subset = region != "Chaco", data = dta); summary(m2) 
se2 <- coeftest(m2, vcov = vcovHC(m2, type = "HC2", cluster = "group")); se2

# Model 3: Prices x Suitability x Settlements, w/o controls
m3 <- plm(ln_conflict ~ I(ln_suitability*ln_price) + I(ln_price*subsistence_farms_median) + I(ln_suitability*ln_price*subsistence_farms_median) + factor(year), index = c("codigo"), model = "within", effect = "individual", subset = region != "Chaco", data = dta); summary(m3) 
se3 <- coeftest(m3, vcov = vcovHC(m3, type = "HC2", cluster = "group")); se3

# Model 4: Prices x Suitability x Settlements, w/ controls
m4 <- plm(ln_conflict ~ I(ln_suitability*ln_price) + I(ln_price*subsistence_farms_median) + I(ln_suitability*ln_price*subsistence_farms_median) + ln_pop + factor(year)*distance_hubs + factor(year)*untitled_share, index = c("codigo"), model = "within", effect = "individual", subset = region != "Chaco", data = dta); summary(m4) 
se4 <- coeftest(m4, vcov = vcovHC(m4, type = "HC2", cluster = "group")); se4

# Model 5: Prices x Suitability x Committees, w/o controls
m5 <- plm(ln_conflict ~ I(ln_suitability*ln_price) + I(ln_price*orgs_locales_median) + I(ln_suitability*ln_price*orgs_locales_median) + factor(year), index = c("codigo"), model = "within", effect = "individual", subset = region != "Chaco", data = dta); summary(m5) 
se5 <- coeftest(m5, vcov = vcovHC(m5, type = "HC2", cluster = "group")); se5

# Model 6: Prices x Suitability x Committees, w/ controls
m6 <- plm(ln_conflict ~ I(ln_suitability*ln_price) + I(ln_price*orgs_locales_median) + I(ln_suitability*ln_price*orgs_locales_median) + ln_pop + factor(year)*distance_hubs + factor(year)*untitled_share, index = c("codigo"), model = "within", effect = "individual", subset = region != "Chaco", data = dta); summary(m6) 
se6 <- coeftest(m6, vcov = vcovHC(m6, type = "HC2", cluster = "group")); se6

# Create table
stargazer(
  m1, m2, m3, m4, m5, m6, 
  se = list(se1[,2], se2[,2], se3[,2], se4[,2], se5[,2], se6[,2]), 
  p = list(se1[,4], se2[,4], se3[,4], se4[,4], se5[,4], se6[,4]),
  omit = c("factor", "Constant", "distance_hubs", "untitled_share", "ln_pop"),
  omit.stat = c("adj.rsq", "ser", "f"),
  covariate.labels = c("Prices$_{t-1}$ $\\times$ Suitability", 
                       "Prices$_{t-1}$ $\\times$ Settlements", 
                       "Prices$_{t-1}$ $\\times$ Suitability $\\times$ Settlements", 
                       "Prices$_{t-1}$ $\\times$ Committees", 
                       "Prices$_{t-1}$ $\\times$ Suitability $\\times$ Committees"),
  add.lines = list(c("","","","","",""), 
                   c("Controls", "No", "Yes", "No", "Yes", "No", "Yes")),
  #c("Municipality FE", rep("Yes", 6)),
  #c("Year FE", rep("Yes", 6))),
  dep.var.labels.include = FALSE,
  model.names = FALSE,
  df = FALSE, 
  multicolumn = TRUE, 
  type = "latex",
  star.cutoffs = c(0.1, 0.05, 0.01),
  star.char = c("*", "**", "***"),
  #float = TRUE, 
  #float.env = "sidewaystable",
  column.sep.width = ".5pt",
  font.size = "scriptsize",
  style = "ajps",  
  title = "Excluding Municipalites from the Chaco Region", 
  label = "chaco",
  notes.append = FALSE,
  notes = c("\\parbox[t]{14.25cm}{\\textit{Note}: All models include municipality and year fixed effects. Standard errors clustered by municipality in parentheses. The sample excludes all the municipalities from the western Chaco forest.}",
            "$^{*}p<.1$, $^{**}p<.05$, $^{***}p<.01$"),
  out = paste0(main,"1_writing/2023/tables/table_wochaco.tex")
)


#### Figure D.2: Marginal Effect of Agricultural Prices on Peasant Resistance by Land Suitability (Excluding the Chaco Region) ----

# Based on model 1
m1 <- plm(ln_conflict ~ ln_suitability*ln_price + factor(year), index = c("codigo"), model = "within", effect = "individual", subset = region != "Chaco", data = dta)

# Covariance matrix for cluster standard errors
covMat1 <- vcovHC(m1, type = "HC2", cluster = "group")

# Plot
pdf(file = paste0(main, "1_writing/2023/figures/inter_wochaco.pdf"), width = 5, height = 5)
interaction_plot_continuous(m1, effect="ln_price", moderator="ln_suitability", interaction="ln_suitability:ln_price", varcov=covMat1, minimum="min", maximum="max", incr="default", num_points = 100, conf=.95, mean=FALSE, median=FALSE, alph=50, rugplot=F, histogram=T, title=NULL, xlabel="Land Suitability", ylabel="Marginal Effect of Agricultural Prices")
dev.off()


#### Figure D.3: Marginal Effect of Agricultural Prices on Peasant Resistance by Land Suitability (Excluding the Chaco Region) ----

# Based on model 3 & 5
m3 <- plm(ln_conflict ~ ln_suitability*ln_price + ln_price*subsistence_farms_median + ln_suitability*ln_price*subsistence_farms_median + factor(year), index = c("codigo"), model = "within", effect = "individual", subset = region != "Chaco", data = dta)
m5 <- plm(ln_conflict ~ ln_suitability*ln_price + ln_price*orgs_locales_median + ln_suitability*ln_price*orgs_locales_median + factor(year), index = c("codigo"), model = "within", effect = "individual", subset = region != "Chaco", data = dta)

# Covariance matrices for cluster standard errors
covMat3 <- vcovHC(m3, type = "HC2", cluster = "group")
covMat5 <- vcovHC(m5, type = "HC2", cluster = "group")

# Plot (a): subsistence settlements
pdf(file = paste0(main, "1_writing/2023/figures/inter_subsist_wochaco.pdf"), width = 5, height = 5)
interaction_plot_threeway_ylim(m3, effect="ln_price", moderator1="ln_suitability", moderator2="subsistence_farms_median", varcov=covMat3, minimum="min", maximum="max", incr="default", num_points = 100, conf=.95, mean=FALSE, median=FALSE, alph=50, rugplot=F, histogram=F, title=NULL, xlabel="Land Suitability", ylabel="Marginal Effect of Agricultural Prices")
dev.off()
# Plot (b): peasant committees
pdf(file = paste0(main, "1_writing/2023/figures/inter_organiz_wochaco.pdf"), width = 5, height = 5)
interaction_plot_threeway_ylim(m5, effect="ln_price", moderator1="ln_suitability", moderator2="orgs_locales_median", varcov=covMat5, minimum="min", maximum="max", incr="default", num_points = 50, conf=.95, mean=FALSE, median=FALSE, alph=50, rugplot=F, histogram=F, title=NULL, xlabel="Land Suitability", ylabel="Marginal Effect of Agricultural Prices")
dev.off()

remove(m1, m2, m3, m4, m5, m6, se1, se2, se3, se4, se5, se6, covMat1, covMat3, covMat5)


#### Table D.3: Peasant Resistance and Distance from Markets ----

# Model 1: Distance to Asunción
m1 <- lm(ln_conflict ~ distance_asuncion + factor(departamento), data = dta_mun); summary(m1)
se1 <- coeftest(m1, vcov = vcovHC(m1, type = "HC0", cluster = "group")); se1

# Model 2: Distance to Export Hubs
m2 <- lm(ln_conflict ~ distance_hubs + factor(departamento), data = dta_mun); summary(m2)
se2 <- coeftest(m2, vcov = vcovHC(m2, type = "HC0", cluster = "group")); se2

# Create table
stargazer(
  m1, m2,
  se = list(se1[,2], se2[,2]), 
  p = list(se1[,4], se2[,4]),  
  covariate.labels = c("Distance from Asunci\\'{o}n",
                       "Distance from Export Hubs"), 
  add.lines = list(c("","","","","")), 
  #c("Department FE", rep("Yes", 5)),
  #c("Year FE", rep("No", 5))),
  omit = c("factor"),
  omit.stat = c("adj.rsq", "ser", "f", "lr", "aic"),
  dep.var.labels.include = FALSE,
  model.names = FALSE,
  df = FALSE, 
  multicolumn = TRUE, 
  star.cutoffs = c(0.1, 0.05, 0.01),
  star.char = c("*", "**", "***"),
  #float = TRUE, 
  #float.env = "sidewaystable",
  font.size = "small",
  #table.placement = "p!",
  style = "ajps", 
  title = "Peasant Resistance and Distance from Markets", 
  label = "distance",
  type = "latex",
  notes.append = FALSE,
  notes = c("\\parbox[t]{8cm}{\\textit{Note}: All models include departmental effects. Huber-White robust standard errors in parentheses. The data come from INE.}",
            "$^{*}p<.1$, $^{**}p<.05$, $^{***}p<.01$"),
  out = paste0(main,"1_writing/2023/tables/table_distance.tex")
)


#### Figure D.4: Local Regressions of Peasant Resistance on Distance to Markets ----

# Plot (a): distance to Asunción
lwr.asuncion <- ggplot(dta_mun, aes(distance_asuncion, ln_conflict)) + 
  geom_smooth(method = "loess", se = TRUE, formula = y ~ x) +
  geom_rug(col = rgb(.5,0,0, alpha=.2), sides = "b") +
  coord_cartesian(xlim = c(0,450)) +
  geom_point() + 
  theme_bw() +
  theme(panel.border = element_blank(),
        axis.text = element_text(size = 14, hjust = 1),
        axis.title = element_text(size = 17, hjust = .5)) +
  xlab("Distance from Asunción") + ylab("Sum of Resistance Events (IHS)")

pdf(file = paste0(main, "1_writing/2023/figures/loess_asuncion.pdf"))
lwr.asuncion
dev.off()

# Plot (b): distance to export hubs
lwr.hubs <- ggplot(dta_mun, aes(distance_hubs, ln_conflict)) + 
  geom_smooth(method = "loess", se = TRUE, formula = y ~ x) +
  geom_rug(col = rgb(.5,0,0, alpha=.2), sides = "b") +
  coord_cartesian(xlim = c(0,450)) +
  geom_point() + 
  theme_bw() +
  theme(panel.border = element_blank(),
        axis.text = element_text(size = 14, hjust = 1),
        axis.title = element_text(size = 17, hjust = .5)) +
  xlab("Distance from Export Hubs") + ylab("Sum of Resistance Events (IHS)")

pdf(file = paste0(main, "1_writing/2023/figures/loess_hubs.pdf"))
lwr.hubs
dev.off()


